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Abstract 

We present results from a comprehensive numerical study of mor- 
phological phase separation (MPS) in unstable thin liquid films on 
a 2-dimensional substrate. We study the quantitative properties of 
the evolution morphology via several experimentally relevant markers, 
e.g., correlation function, structure factor, domain-size and defect-size 
probability distributions, and growth laws. Our results suggest that 
the late-stage morphologies exhibit dynamical scaling, and their evolu- 
tion is self-similar in time. We emphasize the analogies and differences 
between MPS in films and segregation kinetics in unstable binary mix- 
tures. 
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Introduction 



Nonequilibrium processes are of great importance in science and technology. 
The properties of a system are governed by both its constituents and the 
kinetic processes of formation. In this context, there has been intense re- 
search interest in the kinetics of phase transitions, i.e., the nonequilibrium 
evolution of a system which has been rendered thermodynamically unstable 
by a sudden change of parameters, e.g., temperature, pressure, etc. [I]. 

Two important problems in this area are the kinetics of unstable thin 
films [21 HI [5J El [7J and phase- separation kinetics in unstable mixtures 
[HI El QUI EH]- We have recently initiated a study of the analogies and dif- 
ferences between these two physical problems [121 [13] • Both systems are 
characterized by the emergence and growth of domains. In unstable films, 
the coarsening domains consist of flat regions with the equilibrium height h m . 
These are separated by high-curvature regions with gradients in the height 
field. For h m = 0, we have true dewetting (TD) [14J, where holes are punc- 
tured in the film. For h m > 0, the evolution is referred to as morphological 
phase separation (MPS) [151 [TH [17] . O n the other hand, in phase-separating 
binary (AB) mixtures, the growing domains consist of coexisting A-rich and 
i?-rich regions, which are locally in equilibrium. In both systems, coarsening 
is governed by transport of material via diffusive and hydrodynamic pro- 
cesses. The underlying growth mechanisms dictate the system morphology 
and domain growth laws. 

In this communication, we present results from a comprehensive numeri- 
cal study of MPS in 3-dimensional thin films. We use several experimentally 
relevant quantities to completely characterize the emergent morphology. The 
quantities we measure (e.g., correlation function, structure factor, domain- 
size distributions) have often been used to describe evolution morphologies 
in segregating mixtures. However, to the best of our knowledge, these have 
not been used in the context of unstable films. At appropriate places in this 
paper, we will highlight the novel features of MPS in thin films. 

Model and Simulation Details 

The evolution of unstable films is modeled by a continuity equation for the 
height field. This is obtained by simplifying the hydrodynamic equations of 
motion via the lubrication approximation [TSj. In dimensionless units [12J, 
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the continuity equation has the form: 
±H(X,T) = V 



\5H 



(1) 



where H(X,T) denotes the film height at space point X (lying on a d = 2 
substrate) and time T. The unstable initial state of the film is H(X, 0) = 
1+ small fluctuations, i.e., the film is homogeneous to start with. In our 
dimensionless rescaling, all heights are measured in units of the (dimensional) 
initial film height ho- In Eq. (JTJ), the height-dependent mobility M(H) = 
H 3 , corresponding to Stokes flow with no slip. The equilibrium states are 
determined from the free-energy functional: 



F[H] 



dX 



(2) 



where f(H) is the local free energy, and the square-gradient term measures 
the surface tension. We consider a thin film on a coated substrate with 
a long-range van der Waals attraction due to the substrate, and a short- 
range van der Waals repulsion due to the coating [19]. The corresponding 
dimensionless potential is 



f(H) 



1-R 



7 + Jp 



(3) 



6 l(H + D 

In Eq. ([3]), R is the ratio of the effective Hamaker constants, and D is the 
dimensionless coating thickness in units of ho. 

In Fig. [U we plot f(H) vs. H for typical parameter values, R = —0.1 and 
D = 0.2. We also plot f"(H) vs. H in Fig. [TJ The homogeneous thin film 
is spontaneously unstable to fluctuations about the initial state (H Q = 1) 



when f"(l) < 0. Replacing Eqs 
evolution equation: 

dH 



in Eq. (jTJ), we obtain the nonlinear 
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has a wavelength 



3(H 
1 



+ 



R 



D) 3 3H 3 



-V 2 H 



(4) 



shows that the most unstable mode 
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Figure 1: Plot of the free energy f(H) in Eq. ([3]), and its second derivative 
f"(H), with film thickness H. The parameter values are R = —0.1 and 
D = 0.2. The thickness of the equilibrium flat film phase is given by H m = 
J D| J R| 1/3 /[(l + |i?|) 1/3 - \R\ 1/3 ]- For f"(H) < 0, the homogeneous thin film is 
spontaneously unstable and segregates into flat domains and high-curvature 
droplets. 



The final states of the evolution are obtained from the double-tangent 
construction for f(H) in Fig. [TJ The film segregates into phases with H = 
H m = D\R\V*/[(l + \R\) 1/3 - |i?| 1/3 ] and H = oo. The latter value results 
in steepening high-curvature regions which never reach "equilibrium". It 
is relevant to ask whether Eq. (jlj) has a static bump solution H(X), with 
H(X) — > H m as X — >■ ±oo. Such a solution must satisfy the zero-current 
condition: 

g -/'(")-, (6) 

where the constant c = as the bump is flat when X — > ±oo, i.e., d 2 H/dX 2 — > 
0. The first integral of Eq. ()6]) yields 

dH\ 2 



^) =2[f(H)-f(H m )]. (7) 

As H = H m is the sole minimum of f(H) (see Fig. [1]), the right- hand- side 
of Eq. (j7]) can only be zero at H = H m . This immediately rules out a bump 
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solution, which must satisfy dH/dX = at its peak position. The above sce- 
nario should be contrasted with the phase separation of AB mixtures, where 
the system segregates into coexisting A-rich and I?-rich domains. As we dis- 
cuss shortly, this fundamental difference will have important consequences 
for the evolution morphology in thin films. 

We numerically solve Eq. (T4]) in d = 3 (i.e., on a d = 2 substrate), starting 
with a small-amplitude (~ 0.01) random perturbation about the mean film 
thickness H — 1. The parameters D and R were chosen so that the film is 
spinodally unstable at H — 1. The system size is V = (32Lm) 2 - Periodic 
boundary conditions are applied at the lateral ends. A 32-point grid per Lm 
was found to be sufficient when central differencing in space with half-node 
interpolation was combined with Gear's algorithm for time- marching. This 
scheme is especially suitable for stiff equations. 

Numerical Results and Discussion 

In Fig. 121 we show MPS in a thin film evolving from the unstable homo- 
geneous state. The parameter values are R = —0.1 and D = 0.2. The 
snapshots on the left show the height field: regions with H < 1 are un- 
marked, and regions with H > 1 are marked black. The conservation law 
ensures that H = J dXH(X, T) /V is constant in time: H = H = 1 in this 
case. The frames on the right show the variation of the height field along a 
diagonal cross-section of the snapshots [H(X,Y = X,T) vs. X]. The early 
time regime (T = 40) corresponds to the growth of fluctuations about the 
homogeneous state. This growth is exponential and can be obtained from a 
linear stability analysis of Eq. fl3} [20]: H(X,T) = 1 + 6H(X,T), where 

6H(X, T) ~ exp {tV 2 if (1) - V 2 ]} SH(X, 0). (8) 

The growing fluctuations are saturated by the nonlinearity for H ~ H m 
(T = 150): there is no corresponding saturation for the regions with H > 1. 
In the late stages (T = 250), there is growth of domains with H = H m 
(flat phase). At the same time, the defects (or hills) become sharper and 
sharper. Domain growth is driven by the transport of liquid from smaller 
hills to larger hills due to the chemical-potential gradient. This transport 
can be hydrodynamic (when the black regions are connected, e.g., T = 150) 
or diffusive (when the black regions are not connected, e.g., T = 250). 
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Figure 2: Kinetics of morphological phase separation (MPS) in an un- 
stable thin film in d = 3. The initial condition for the evolution con- 
sisted of uniformly-distributed small-amplitude random fluctuations about 
H(X,Y,T = 0) = 1. The parameter values were R = —0.1 and D = 0.2. 
Other simulation details are provided in the text. The dimensionless evo- 
lution times are specified above each frame. Left: Snapshots of the height 
field. Regions with H > 1 (high-curvature phase) are marked in black, 
whereas those with H < 1 (flat phase) are unmarked. Right: Variation of 
the height field along a diagonal cross-section (Y = X) of the snapshots. 



Notice that the morphology evolves continuously in Fig. |2j from bicon- 
tinuous (T = 40) to circular domains of flat phase (T = 150) to droplets of 
high- curvature phase (T = 250). We expect that the system enters a scaling 
regime for T > 250, where the morphology becomes self-similar and only the 
domain size grows. A similar evolution is seen for much thicker films with 
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D = 0.8 in Fig. [3J As expected, the time-scales of MPS are much longer for 
the thick film. The statistical results we show subsequently correspond to 
the film with D = 0.2, whose evolution is shown in Fig. |2j 





Figure 3: Analogous to Fig. [2} but for the parameter values R = —0.1 and 
D = 0.8 (corresponding to a much thicker film). 

Let us present some simple arguments to understand the evolution of the 
scaling morphology in Fig. [2] (e.g., T = 250) and Fig. [3] (e.g., T = 9.7 x 10 5 ). 
Let L(T) and Hd(T) denote the characteristic size and height of a droplet 
of the high- curvature phase at time T. If the number of droplets in the 
system is N(T), the amount of surplus liquid in the defects ~ NL 2 Hd- The 
corresponding depletion of liquid in the flat domains is (V — NL 2 )(1 — H m ), 
where V is the substrate area. As there is a conservation of liquid, we have 

NL 2 H d ~(V-NL 2 )(l-H m ), (9) 
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or 

6(T) ~ (1 ~ Hm) ~ (10) 

In Eq. ffTOj) . 0(T) = N(T)/V denotes the number of droplets per unit area. 
In the final step of Eq. fflQj) . we have neglected (1 — H m ) (which is constant) 
with respect to H d (which grows). 

How do these droplets grow with time? Their coarsening proceeds via an 
evaporation-condensation mechanism. The smaller droplets evaporate and 
their material is diffusively transported to the larger droplets through the 
flat domains. The chemical potential at a point is estimated as 

^ = f{H)-V 2 H. (11) 

We assume that the droplets have a conical profile (cf. Fig. [2]) with the height 

H{r) ~ H d - (H d - 1) ( z ) , (12) 

where the droplet is centered at r = 0, and n is the profile exponent. Then 

v^Ii^W^^! (13) 

r or \ or J rL 

Thus, we estimate the chemical potential at the droplet boundary (H = 1) 

as 

^ m + (JkzJK^ m + Ef. (14) 

A better estimate of the chemical potential can be obtained by integrating 
over the droplet. The resultant value is the same as that in Eq. (TH|) up to 
geometric prefactors. 

The chemical-potential gradient occurs over flat domains (with H = H m ) 
of size ~ L, and the corresponding current is 

J ~ -M(H m )\Vfi\ ~ . (15) 

Thus, the droplets grow as 

alL HlH d n 2 

dT ~ I? ' (16) 
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In the scaling regime, we expect Hd ~ L so that L(T) ~ (HlT) L/ . This 
growth law is well known in the context of phase-separation kinetics, and is 
referred to as the Lifshitz-Slyozov (LS) growth law [H [21] . In earlier work 
[T2"| [T5] , we have demonstrated that MPS in the d = 2 film (i.e., on a d = 1 
substrate) also obeys the LS growth law. The above arguments apply for 
d = 2 films also, except when n — 1 in Eq. ( 112)) . 

There are several experimental tools to quantitatively characterize the 
MPS evolution morphologies. The correlation function of the height field is 
defined as 

C(r,T)=C(r,T) = (SH(X,T)SH(X + f,T)>, (17) 

where 5H (= H — i/ ) denotes the fluctuation in the height field. The 
angular brackets in Eq. ( !T7|1 denote an averaging over independent runs. 
The statistical results presented here are obtained as an average over 10 runs. 
The correlation function depends only on the distance r (the magnitude of 
r) because the system is translationally invariant and isotropic. In scattering 
experiments with, e.g., light, X-rays, neutrons, etc., we measure the structure 
factor, which is the Fourier transform of C(f, T) with wave- vector A;: 

S(k,T) = S{k,T) = J dre tlr V(r,T). (18) 

Let us present numerical results for the statistical quantities defined 
above. In Fig. HJ^a), we plot the spherically-averaged correlation function 
[C(r,T)/C(0,T) vs. r] for the evolution depicted in Fig. [2j At early times, 
linear theory applies and the corresponding structure factor is [2D] 



Siin(Jfe, T) ~ AV exp [2k 2 (a - k 2 ) T] , (19) 

where A is the amplitude of initial fluctuations, (SH (X, 0)5H(X', 0)} = 
A5(X— X'). In Eq. f[T9"|) . V denotes the system volume, and a = —f"{l). The 
early-time correlation function is obtained as the inverse Fourier transform 
of S\ in (k,T): 

CW(r,T) = AV J ^Le-^exp [2k 2 (a - k 2 ) T] 
AV f°° 

= dk k.J (kr) exp [2k 2 (a - k 2 ) T] . (20) 

2vr Jo 

The solid line in Fig. Ufa) denotes the expression for Cn Q (r,T)/Cu n (Q,T) 
from Eq. ( 12"U|) with T = 40. It is in excellent agreement with the numerical 
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data for T = 40. As time goes on, the fluctuations grow exponentially and 
linear theory ceases to hold. There is no change in the length scale during 
the exponential growth regime - hence, the zero-crossings of the correlation 
function at T = 150 in Fig. 11(a) are comparable to those for T = 40. (The 
length-scale data shown later demonstrates that domain growth occurs for 
T > 150.) There is a crossover in C(r, T), as expected from the morphological 
evolution seen in Fig. [21 The snapshot at T = 250 is indicative of the scaling 
morphology, viz., droplets of the high-curvature phase in a background of 
the flat phase. In the late stages, we expect dynamical scaling of C(r, T) and 



where L is the characteristic length scale. In Eq. ( 121]) . g(x) and f(p) are 
scaling functions which do not depend on time. In Fig. H^a), the correlation 
function at T = 250 lies in the scaling regime. 

In Fig.Hfb), we plot the spherically-averaged structure factor [S(k,T) vs. 
k, on a log-log scale] for the same times as in Fig. 0Ja). As before, there 
is a crossover in the functional form of the structure factor. The data set 
for T = 40 corresponds to the linear regime, and is in good agreement with 
the expression in Eq. (fl9|) . Notice that the linear theory is not valid at very 
large values of k (> 0.5) where we see effects of discreteness of the simulation 
lattice. The data set at T = 250 corresponds to the asymptotic scaling form 
in Eq. (T2TT) . The peak location [k m (T) ~ moves to smaller values as the 
characteristic length scale grows. There is a shoulder at higher values of k, 
which reflects the second length scale apparent in the snapshots in Fig. [2j 
This shoulder disappears as we approach the asymptotic droplet morphol- 
ogy. Finally, the tail of the structure factor does not show the Porod tail 
[23], S(k,T) ~ ^-(d+i) a t i ar g e ^ which is characteristic of phase-separating 
systems |21j. The Porod tail arises from scattering off sharp interfaces - how- 
ever, the MPS morphology in Fig. [2] does not have any equilibrium interfaces 
at all. 

An alternative method of characterizing the growth morphology is via the 
probability distribution P(l, T) of domain size I. Consider a typical snapshot 
in Fig. [21 We examine the variation of the height field in both the x- and y- 
directions, and keep track of the "zero" -crossings (where H — 1, the average 
height). This yields the distribution of sizes for the flat domains, as well as the 



S(k,T) [221 H]: 



C(r,T) 
S(k,T) 



9(r/L), 
L d f(kL), 



(21) 
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Figure 4: Statistical quantities for the evolution depicted in Fig. |2j These 
are obtained on (32Lm) 2 systems as an average over 10 independent runs, 
(a) Plot of the spherically-averaged correlation function C(r,T)/C(0,T) vs. 
r. The solid line denotes the result from linear theory in Eq. ( 120]) . evaluated 
at T = 40. (b) Plot of the spherically-averaged structure factor S(k, T) vs. 
k. The solid line denotes the expression in Eq. f fl9|) . evaluated at T = 40. 
The line with slope —3 refers to Porod's law in d = 2, S(k, T) ~ fc-id+i) f or 
large k. The Porod law characterizes scattering off sharp interfaces [23l [21] . 
In the unstable thin film, there are no equilibrium interfaces as there is only 
one homogeneous phase. 



defects or hills. In Fig. [5(a), we plot the domain-size distribution Pdam{l,T) 
vs. I. This quantity also shows the expected crossover from early to late 
times. The data set for T = 40 has a broader distribution than the data set 
for T = 150, though the peak positions are similar. The T = 40 morphology 
has a broader distribution because it still contains modes with various length 
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scales - these are dominated by the most unstable wavelength (L M ) as the 
fluctuations grow. The peak positions at T = 40, 150 are comparable because 
there is no domain growth in the early exponential regime. At later times, 
the domains coarsen and the peak position of Pdo m (l,T) shifts to the right 
[see data set for T = 250 in Fig. |5]^a)]. The functional form at T = 250 
corresponds to the asymptotic scaling function Pdom(a ; ), which is defined as 

P dom (/,T) = L~ l p dora (l/L). (22) 

We make two further observations about P dom (l,T). First, the plots at 
early times (T = 40, 150) show a marked shoulder, consistent with our obser- 
vation of a two-scale morphology in the snapshots of Fig. |2j and the structure 
factors of Fig. BJb). This shoulder disappears as we enter the asymptotic 
scaling regime. Second, the tail of the distribution decays exponentially, i.e., 
linearly on the semi-log plot in Fig. [5]^a). This is typical of morphologies 
with well-defined characteristic scales, and has been observed earlier in the 
context of phase-separation morphologies [24] . 

In Fig. Mjo), we plot the defect-size distribution [P de f(l,T) vs. I] at the 
same times as in Fig. |5]^a). This plot shows the same general features as 
Fig. da). 

Finally, in Fig. |H1 we study the growth of the characteristic length scale. 
We use three different measures of the length scale, obtained from the sta- 
tistical quantities introduced earlier: 

(a) Lc, the length scale up to which C(r,T) decays to half its maximum 
value (which arises at r = 0). 

(b) Ls, the inverse of the first moment of S(k,T): L$ = (k)^ 1 with 

_ J-dkkS(k,T) 
{k) ~ J °°dkS(k,T) • {26) 

(c) Lp, the average domain size from the relevant probability distribution: 

POO 

L P = / dllP dom (l,T). (24) 
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We show the time- dependence of these length scales in Fig. [6j It is interesting 
to note that all the different length scales show the same physical behavior. 
There is an early regime (up to T ~ 150), where the initial fluctuations grow 
exponentially but there is almost no change in the length scale. In the context 
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50 100 150 

I 

Figure 5: Probability distributions of domain and defect sizes for the evolu- 
tion depicted in Fig. |2j The data is plotted on a linear-log scale. The solid 
lines denote exponential fits to the tail region, (a) Domain-size distribution, 
-Pdom(^,^) vs. I. (b) Defect-size distribution, P^{l,T) vs. I. 

of phase-separation kinetics, this is referred to as the Cahn-Hilliard-Cook 
(CHC) regime. This is followed by an intermediate regime, where nonlinear 
effects saturate the growing fluctuations. The system segregates into domains 
of the flat phase and defects of the high-curvature phase. In the late stages 
(T > 175), there is growth of domains. It is not possible to ascertain the 
asymptotic growth exponent from our present results, which only access a 
small window in the scaling regime. For this, we require 1-2 further decades 
of growth, which would be computationally very expensive. The numerical 
difficulty in accessing the late stages can be understood as follows. Recall that 
the interfaces become progressively steeper with time. Therefore, we need a 
corresponding reduction of the mesh size in space (and time) to resolve the 
interfaces properly. 
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Figure 6: Time-dependence of the characteristic length scale for the evolution 
in Fig. [2J We show data for three measures of the length scale: Lc, L$, and 
Lp, obtained from C(r,T),S(k,T), and Pdom{l,T), respectively. 



Summary 

Let us conclude this paper with a brief summary. We have undertaken a 
comprehensive numerical study of morphological phase separation (MPS) in 
unstable liquid films. In particular, we have used several experimentally rele- 
vant tools to quantitatively characterize the emergent morphologies. The sta- 
tistical quantities we measure are the correlation function, structure factor, 
domain-size and defect-size probability distributions, and the corresponding 
growth laws. The study of these properties provides a complete picture of the 
evolution dynamics. We make several important predictions in this context. 
At late times, the system should enter a scaling regime, where the above 
quantities show dynamical scaling, i.e., the morphology becomes self-similar 
in time. We hope that our numerical study will motivate fresh experiments 
on unstable thin films. 

There are important analogies and differences between MPS in thin films 
and phase-separation kinetics in binary mixtures. We have highlighted some 
of these features in this paper. A crucial difference between the two systems 
is the nature of the coarsening domains and defects. The unstable film segre- 
gates into droplets of a high-curvature phase (defects), separated by domains 
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of a flat phase. On the other hand, the unstable AB mixture segregates into 
coexisting domains of A-rich and I?-rich phases, which are separated by near- 
equilibrium interface defects. This difference has important consequences for 
the evolution morphologies and their quantitative properties. 

R.K. acknowledges the financial support of the Department of Science 
and Technology, India. 
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